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ABSTRACT 

Aims. In this paper, we explore the diagnostic power of the far-IR fine-structure lines of [Oi] 63.2/jm, 145.5 /jm, [Cn] 157.7 //m, as 
well as the radio and sub-mm lines of CO J=l-0, 2-1 and 3-2 in application to disks around Herbig Ae stars. We aim at understanding 
where the lines originate from, how the line formation process is affected by density, temperature and chemical abundance in the 
disk, and to what extent non-LTE effects are important. The ultimate aim is to provide a robust way to determine the gas mass of 
protoplanetary disks from line observations. 

Methods. We use the recently developed disk code ProDiMo to calculate the physico-chemical structure of protoplanetary disks and 
apply the Monte-Carlo line radiative transfer code Ratran to predict observable line profiles and fluxes. We consider a series of Herbig 
Ae type disk models ranging from 10~ 6 M Q to 2.2x 10~ 2 M Q (between 0.5 and 700 AU) to discuss the dependency of the line fluxes 
and ratios on disk mass for otherwise fixed disk parameters. This paper prepares for a more thorough multi-parameter analysis related 
to the Herschel open time key program Gasps. 

Results. We find the [Cn] 157.7/mi line to originate in LTE from the surface layers of the disk, where T„ # 7j. The total emission 
is dominated by surface area and hence depends strongly on disk outer radius. The [Oi] lines can be very bright (> 10~ 16 W/m 2 ) and 
form in slightly deeper and closer regions under non-LTE conditions. For low-mass models, the [Oi] lines come preferentially from 
the central regions of the disk, and the peak separation widens. The high-excitation [Oi] 145.5yum line, which has a larger critical 
density, decreases more rapidly with disk mass than the 63.2 fim line. Therefore, the [Oi] 63.2 yum/145.5 fim ratio is a promising disk 
mass indicator, especially as it is independent of disk outer radius for R out > 200 AU. CO is abundant only in deeper layers A v <: 0.05. 
For too low disk masses (M disk S 10~ 4 M Q ) the dust starts to become transparent, and CO is almost completely photo-dissociated. For 
masses larger than that the lines are an excellent independent tracer of disk outer radius and can break the outer radius degeneracy in 
the [Oi] 63.2 fjm/[C n] 157.7 fan line ratio. 

Conclusions. The far-IR fine-structure lines of [Cn] and [Oi] observable with Herschel provide a promising tool to measure the disk 
gas mass, although they are mainly generated in the atomic surface layers. In spatially unresolved observations, none of these lines 
' carry much information about the inner, possibly hot regions < 30 AU. 
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'jjj . 1. Introduction from the dust spectral energy distribution (S EP) are generally 

smaller than those derived from the gas l ine dlsella et all I2007I : 
Observations of gas in protoplanetary disks are intrinsically dif- |Pietu et all I2005I) . iHughes et al] d2008al) suggest a soft outer 
ficult to interpret as they reflect the interplay between a complex e d ge as a solution to the discrepancy. Comparison of CO J=3- 
chemical and thermal disk structure, statistical equilibrium and 2 maps of four disks to different types of dis k models strongly 
optical depth effects. This is particularly true if non-thermal ex- supports a soft edge in favor of a sharp cutoff. lPietu et alj f2007) 
citation such as fluorescence or photodissociation dominate the use the CO and HCO + lines to probe the radial and vertical tern- 
statistical equilibrium. perature profile of the disk. Simple power law disk models and 
The first studies of gas in protoplanetary disks concentrated LTE radiative transfer provides best matching results for radial 
on the rotational transiti ons of abundant molecules such a s temperature gradients around r~ 0t5 . The 12 CO J=2-l line, 13 CO 

J=2-l and 13 CO J=l-0 lines used in their analysis probe subse- 



CO, HCN and HCO + (e.g . 
7993t iDutteTetail 1997: 



Beckwith et al., 1986; Koerner et al 



van Zadelhoff et all I200U iThi et al 



quently deeper layers and reveal a vertical temperature gradient 

20041) . Those lines originate in the outer regions of disks, r > ranging from 50 K in the higher layers to below the freeze-out 

100 AU, where densities are at most n ~ 10 7 cm 3 . The inter- temper ature of CO in the m idplane. This confirms earlier find- 

pretation of those lines was mainly based on tools and exper- ings bv lDartois et all d2003l) . 

tise developed f or molecular clouds. Using the CO J=3-2 line, However, disk masses derived from the CO lines are in gen- 

iDent et all d2005h inferred for a sample of Herbig Ae and Vega- eral lower than disk ma s ses derived from dust observations (e.g. 

type stars a trend of disk outer radius with age; on average, the IZuckerman et all 119951: IThi et all 1200 ll) . Possible explanations 

outer disk radius in the 7-20 Myr range is three times smaller include CO ice formation in the cold midplane and photodisso- 

than that in the < 7 Myr range. Also, the disk radii inferred ciation in the upper tenuous disk layers. Any single gas tracer 
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alone can only provide gas masses of the species and volume 
from which it originates, the same way as dust masses derived 
from a single photometric measurement are only sensitive to 
grains of a particular size range, namely those grain sizes that 
dominate the emission at that photometric wavelength. Hence 
individual gas tracers are more valuable for probing the physi- 
cal conditions of the volume where they arise than the total disk 
mass. A combination of suitable gas tracers can then allow us to 
characterize the gas properties in protoplanetary disks and study 
it during the planet formation process. 

This paper aims at exploring the diagnostic power of the fine 
structure lines of [O i] and [C n] in the framework of upcoming 
Herschel observations. Earlier modeling of these lines indicated 
that they should be detectable down to disk masse s of 10~ 5 Mr 
of ga s, so also in the very gas-poor debris disks dKamp et al.L 
2005). More recent work by lMeijerink et al.l (120081) presents fine 
structure lines from the inner 40 AU disk of an X-ray irradiated 
T Tauri disk; the models indicate that the [Oi] emission origi- 
nates over a wide range of radii and depth and is sensitive to the 
X-ray luminosity. However, most of the C n and O i line emis- 
sion comes from larger radii wh ere the gas temperatur e is domi- 
nated by UV heating processes. iJonkheid et al.1 d2007) find from 
thermo-chemical models of UV dominated Herbig Ae disks that 
the [Oi] lines are generally a factor 10 stronger that the [Cn] 
lin e. We started to expl ore the origin of the fine structure lines 
in IWoitke etall d2009al) and find that the [Cn] 157.7 fim line 
probes the upper flared surface layers of the outer disk while the 
[Oi] 63.2 yum line originates from the thermally decoupled sur- 
face layers inward of about 100 AU, above Ay ~ 0.1. The latter 
line is very sensitive to the gas temperature and might be used to 
distinguish between hot (r gas * 1000 K) and cold (r gas = T^st) 
disk atmospheres. Since the fine structure lines generally orig- 
inate from a wider radial and vertical range than for example 
the I2 CO rotational lines, they are potentially more suitable gas 
mass tracers. 

We use in this paper the di sk modeling code ProDiMo 
presented in Woit ke et al.l (l2009al) to study the gas line emis- 
sion from disks around Herbig Ae stars. The main focus are 
the fine-structure lines of Cn and Oi which will be observed 
for a large sample of disks during the Herschel open time 
Key Program Gasps (Gas evolution in protoplanetary systems: 
http://www.laeff.inta.es/projects/herschel). The disk parameters 
were chosen to resem ble the disk around MWC480. Ac c ording 
to p revious wo r k bvlMannings et al.1 dl997l) . iThi et al.1 d2001l) 
and iPietu et ail d2007l) . the disk around this star extends from 
« 0.5 AU to 700 AU. We choose here a surface density profile 
E ~ r~ . The central star is an A2e Herbig star with a mass of 
2.2 M and an effective temperature of 8500 K. 

Section [2] gives a short summary of the disk modeling ap- 
proach. The line radiative transfer method, re-gridding and the 
atomic input data are described in Sect. [3] We then briefly 
discuss some basic properties of the Herbig Ae disk models 
(Sect. HJ before we present the fine-structure lines (Sect. |5]l 
and conclude with a discussion of the diagnostic strength of 
fine structure line ratios and a comparison to previous ISO and 
submm observations (Sect. [6J. 

2. ProDiMo 

ProDiMo is a recently developed code for computing the hy- 
drostatic structure of protoplanetary disks. This code combines 
frequency-dependent 2D dust continuum radiative transfer, ki- 
netic gas-phase and UV photo-chemistry, ice formation, and 
detailed non-LTE heating & cooling with the consistent calcu- 



lat ion of the hydrostat ic disk structure. Details can be found 
in IWoitke et al.l d2009al) . We summarize in the following some 
other aspects that are particularly rel evant to this study. 

We use a Phoenix stellar model iBrott & Hauschildl d2005l) 
with an effective temperature of 8500 K for the stellar irradia- 
tion and a highly diluted 20000 K black body for the IS radia- 
tion field that penetrates the disk from all sides. The dust opac- 
ities are computed u sing Mie theory and optical constants from 
iDraine & Led d 19841) . The chemical network contains 71 species 
(build from 9 elements) connected through 950 reactions (photo 
reactions, CR ionization, neutral-neutral, ion-molecule, as well 
as grain adsorption and desorption processes for CO, CO2, H2O, 
CH4 and NH3 ice). The gas temperature follows from an ex- 
tensive heating and cooling balance that includes fine-structure 
line cooling as well as molecular and optical lines (fluorescence 
in the inner disk). The code does not yet include X-ray heat- 
ing. This seems a minor issue for Herbig stars that generally 
have quite moderat e X-ray luminosities, Lx <s 10 30 erg s _1 (e.g. 
IStelzeret al.ll2006l) . 

The treatment of photoionization and photodissociation in 
this work differs from that in the original code. The photoion- 
ization and photodissociation rates, R ph , are now computed at 
each grid point using the spectral photon energy density Aux cal- 
culated by the 2D continuum radiative tra nsfer and the cross- 
sectio ns cr(X) from the Leiden database dvan D ishoeck et all 
2008). The photorate for continuous absorption is 

R ph = i j cr(A)Au A dA . (1) 

When the photoprocess is initiated by line absorption, the rate 
becomes 

2 

* Ph = ^a^fm^lh) , (2) 

where fj is the oscillator strength for absorption from lower level 
i to upper level j, rjj is the efficiency of state j with values lying 
between and 1, and ne 2 /mc 2 is 8.85 x 10~ 21 with A in A. The 
2D UV radiative transfer uses now three bands defined between 
91.2, 111, 145, and 205 nm with central wavelengths of 100, 
127, and 172 nm and the radiation field at intermediate wave- 
lengths is recovered through spline interpolation from the three 
bands. 

We compute a series of Herbig Ae disk models with masses 
between 2.2 10~ 2 and 10~ 6 M . The parameters are summa- 
rized in Table Q] The dust used in this model is typically larger 
than ISM dust and generates the ~ A~ l opacity law as de- 
rived from dust observations (e.g. [Beckwith & Sargentl ] 199 1 ; 
iMannings & Emersorlll994lRodmann et al.Ll2006l) . 

3. Line radiative transfer 

3.1. Methods 

We use the two-dimen sional Monte Carlo radiative t ransfe r code 
Ratran developed by Hogerheij de & van der Takl d2000h . The 
code uses a two-step approach to solve the non-LTE line ra- 
diative transfer, Amc, and Sky. The first code solves the level 
population numbers for a given model atom/molecule within an 
arbitrary two-dimensional density and temperature distribution. 
The second one performs the ray tracing to derive the emission 
for a given line, distance and disk inclination. In the following, 
we discuss some aspects that are particularly relevant in apply- 
ing these codes to complex chemical disk stratifications. 
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Table 1. Parameters of the Herbig Ae model series. 



Quantity 


Symbol 


Value 


stellar mass 




2.2 M G 


effective temperature 


T c ft 


8500 K 


stellar luminosity 


u 


32 L 


disk mass 


A/disk 


2.2 x 10~ z , 10^, 






10" 3 , 10- 4 , 10" s , 






10- 6 M G 


inner disk radius 




0.5 AU (1) 


outer disk radius 


"out 


700 AU 


radial column density power index 


6 


1.0 


dust-to-gas mass ratio 


PdlP 


0.01 


minimum dust particle radius 


^min 


0.05 yum 


maximum dust particle radius 




90(1 urn 


dust size distribution power index 


^pow 


3.5 


dust material mass density 


Pgr 


2.5 gcirT 3 


strength of incident ISM UV 


* 1SM 


1 


cosmic ray ionization rate of H2 


£cit 


5 x 10~ 17 s~ ] 


abundance of PAHs relative to ISM 


/pah 


0.1 


a viscosity parameter 


a 


0.0 



(1): soft inner edge applied, see Sect. 3.1 o f lWoifkeetal1 ll2009a) 



3.2. Re-gridding 

The 2D non-LTE line transfer code Ratran requires a grid of 
rectangular cells in cylindrical coordinates r e {r,-, r, + i) and 
Z e \Zj,Zj+\) which is different from the grid of points used in 
ProDiMo. Therefore, we have to create a suitable grid of cells 
for Ratran and "fill" the cells in a physically sensible way, which 
will involve some kind of averaging for the physical quantities. 

We choose to evaluate these mean values by integration and 
define the following general function between the ProDiMo grid 
points 



f(r,z) = f rPexp(-(z/H) 2 ). 



(3) 



which can be applied to any physical quantity, e. g. the species 
particle density n sp (r, z), or the product of hydrogen nuclei den- 
sity and gas temperature (n sp ■ r gas )(r, z). Given any point (r,z) 
inside the ProDiMo grid, we determine the indices of the sur- 
rounding 2x2 ProDiMo corner grid points, and fix the free co- 
efficients /o, p and H to fit the values at the corner points. 

Next, we calculate the following integrals over the Ratran 
cells 



Vij = n(rf +1 - rf)(z j+i - Zj ) 

jV sPj ij = In J r J n s? (r,z)dzdr 

n zj 
n+\ z .i+i 

(N S pT m )ij = 2n J r J(n sp -T m )(r,z)dzdr 

n Zj 

(N(H)T iust )ij = 2n J r j '(n (H > ■ r dust )(r, z) dzdr 



(4) 
(5) 

(6) 

(7) 



to derive mean values for the species density, gas and dust tem- 
perature in our Ratran cells 



Tgas jj — vVspTgas/i.7 / ^sp, ij 

7dust,y = (N(H)Tdust)ij / N(u),ij 



(8) 
(9) 
(10) 




r[AU] 

Fig. 1. Distribution of cells for the Ratran radiative transfer 
code. The background contours show the 1 density distribution 
of the 10~ 2 M model. The black boxes indicate the distribution 
of 400 cells across the model (original grid points are denoted 
by black dots). 



Similar formulae apply to the collision partner densities. We 
have choosen this procedure (i) to assure total and species mass 
conservation, (ii) to conserve the total line emission in optically 
thin LTE at long wavelengths (Rayleigh-Jeans approximation) 
and (iii) to guaranty that the total thermal dust emission in the 
optically thin case (Rayleigh-Jeans approximation) is conserved. 

Fig. [TJ shows how the Ratran grid compares to the original 
grid and underlying Oi density distribution (10~ 2 M model). 

3.3. Modified radiative transfer code Ratran 

The original Ratran code uses an accelerated Lambda 
iteration scheme to acceler ate the convergence 
( Hogerheii de & van der Tall 120001) . For rather complex 
molecules with a large number of levels heavily intercon- 
nected through lines of various strengths, such as water, the 
performance is too slow to allow a full exploration of the disk 
parameter space. We describe in the following the implemen- 
tation of two additional methods that significantly improve the 
Ratran performance. 

The line transfer simulations in Ratran are carried out in two 
stages, called the "fixset" and the "random" stages. During the 
fixset stage the number, starting points and directions of the rays 
are fixed, and only the non-local feedback between level popu- 
lations and spectral intensities are solved iteratively. In order to 
accelerate the convergence during the r andom-noise -free fixset 
stage, we have included the procedure o flAueti 0984) for strictly 
convergent transformations, known as the Ng-iteration. During 
the random phase, the number of rays is successively increased 
in all cells until three different sets of rays give approximately 
the same results. 

The typical errors of the Monte-Carlo method depend on 
the chosen set of random numbers. For standard (pseudo- 
random) number generators, the er rors decrease with t he num- 
ber of photon packages as AT 1 / 2 (iNiederreitett 119921) . Ouasi- 
random numbers or low discrepancy sequences, to which the 
Sobol sequence belongs, try to samp l e more uniformly t he dis - 
tributed points dSobol & Shukhmanl l2007t iPress et all 120021) . 
The asymptotic errors of this quasi Monte-Carlo simulations de- 
crease as (log(N)) d /N, where d is the dimension of the prob- 
lem. For Ratran, the dimension is d = 5, because the posi- 
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tion in the cell, the direction and the Doppler velocity shift of 
the photon package are chosen randomly. The _useof quasi- 
random sequences was already advocated by lJuvelal (1 19971) . 
since the gain in number of photon packages is significant, 
but has not been widely adopted. One of the main critics of 
quasi-random number se quences is th at there are only a lim- 
ited number of sequences (ISloanl ll993l). Hybrid methods called 
randomized or scrambled quasi-random sequences combine the 
advantages of both individual methods, infinite number of se- 
quences and low error, resulting in an asymptotic noise that 
decreases in N~ s (los(N)) ^ d ~ 1 ^ 2 , where s is between 1 and 1.5 
dHickernell & Yuel [2000). Scramble quasi-random sequences 
show the fastest noise decrease among the three types of se- 
quences. 

We replaced the pseudo-random sequence in Ratran by the 
quasi-random Sobol sequence combined with a Owen-Faure- 
Tezuka type of scrambling dOwenl 120031) . A detailed descrip- 
tio n of the ACM algori thm 659 implemented in Ratran is given 
by lBratley et al.l dl994l) . Difficult problems, like those involving 
water lines, require up to 10 8 photon packages with the stan- 
dard random sequence but only between 10 4 (s = 1.5) and 
~ 2 x 10 6 {s = 1) photon packages with the scrambled quasi- 
random sequence, a gain in CPU time between 50 (s — 1) and 
10 4 (s = 1.5). This speed gain opens the possibility to study 
water line emission s from protoplaneta ry disks with reasonable 
computing times dWoifke et alll2009bl) . 

We extended the original code also to include a larger num- 
ber of collision partners. This proofs to be important as in the 
original code, the density of the first collision partner is used to 
compute the dust continuum emission. However, the density of 
one of the collision partners does not always equal the total hy- 
drogen number density, thus causing inconsistencies either in the 
dust emission or the collision rates. Moreover, some lines orig- 
inate over a wide range of physical and chemical conditions, so 
that the dominant collision partner changes between the various 
disk regions (e.g. from atomic hydrogen and electrons in the low 
density regions to molecular hydrogen and electrons deeper into 
the disk). 



3.4. Ray tracing 

We assume for our generic disk a distance of 131 pc and and 
inclination of 45 degrees. The pixel size (spatial resolution) is 
0.05" and the entire box has a size of 13 X 13" to ensure that no 
emission is lost (disk diameter is 10"). The velocity resolution 
is set to 0.05 km/s and the total range is from -25 to 25 km/s 
for the C n line and -40 to 40 km/s for the O i lines. The different 
velocity ranges reflect the differences in radial origin of the lines. 
For the oxygen fine structure lines, oversampling of the central 
pixels (out to 1.6") was used, i.e. an additional 2 rays generated 
per pixel. 

The CO lines were computed within the same box of 13 x 
13", but with a spatial resolution of 0.12" and a spectral resolu- 
tion of 0.2 km/s. No oversampling was used in this case. 



3.5. Atomic and molecular data 

Energy levels, statistical weights, Einstein A coefficients and 
collision cross sections are taken from the Leiden Lambda 
database (Scho ier et all [2005). Table |2]provides an overview of 
the C II, O I and CO data used in this work. 



Table 2. Atomic and molecu lar data taken from the Lambda 
database dSchoier et al .1120051) . 

Species # Lev. # Lines Collision Reference 
C + 2 1 = 



3 3 



CO 26 25 



partner 


H 


Launay & Roueff (1977) 


e~ 


Wilson & Bell (2002) 


o-H 2 


Flower & Launay a (1977) 


P-H 2 


Flower & Launay (1977) 


H 


Launay & Roueff (1977) 


e" 


Belletal. (1998) 


o-H 2 


Jaquet et al. (1992) 


p-H 2 


Jaouet etal. (1992) 


H + 


Chambaud (1980) (1) 


H 


Chu & Dalearno (1975) 


H 2 


Schinkeetal. (1985) 



(1): private communication 



Table 3. Atomic and molecular data for the oxygen and carbon 
fine structure lines and the CO rotational lines. 





Line identification 


E„ 


Ei 


A*, 


«crit 






[/an] 




[K] 


[K] 


[s- 1 ] 


[cm- 3 ] 




Oi 


63.2 


J P,- 3 P, 


227.7 


0.0 


8.865(-5) 


5(5) 




Oi 


44.1 


3 Po- 3 P 2 


326.6 


0.0 


1.275(-10) 


0.5 




Oi 


145.5 


3 Po- 3 P. 


326.6 


227.7 


1.772(-5) 


6(4) 




Cn 


157.7 


2 P, - 2 P, 


91.2 


0.0 


2.300(-6) 


3(3) 






Line identification 


E u 


Ei 


A„i 


"ci-it 






[GHz] 




[K] 


[K] 


[s- 1 ] 


[cm" 3 ] 




CO 


115. 3 


J = 1 - 


5.53 


0.0 


7.203(-8) 


« 5.0(3) r- 


-0.66 


CO 


230.5 


J = 2 1 


16.60 


5.53 


6.910(-7) 


(1) 1.9(4) T g 


-0.45 


CO 


345.8 


J = 3-2 


33.19 


16.60 


2.497(-6) 


m 4.6(4) f E 


-0.35 



(1): CO critical densities from lKamp & van Zadelhoff (2001) 
The notation 5(5) stands for 5 10 5 . 

Oxygen levels 2, 1, correspond to 3 P 2 , 3 P[, and 3 Po, respectively. 



3.6. Dust opacities 

The dust opacities used in the radiative transfer code are con- 
sistent with the opacities used in the computation of the disk 
model with ProDiMo. The choice of opacities impacts the con- 
tinuum around the line, i.e. the dust thermal emission. For oxy- 
gen, the low fine structure levels can be pumped by the ther- 
mal dust background. We see differences of up to 50 % in the 
continuum fluxes around the O i fine structure lines and up to 
20 % differences in the line emission itself when we choose 
either the grain size distribut i on fro m Table [1] using optical 
constants from | Draine & Led ([1984) or opacity tables from 
lOssenkopf & Henningl dl994l) . The latter are based on an MRN 
size distribution for interstellar medium grains f(d) ~ a~ 35 with 
sizes 5 nm < a < 0.25 /vm. In this paper, the possibility of icy 
grain mantles or non-spherical shapes is not explored. 

4. The disk models 

We compute a series of six disk models with different disk 
masses (2.2 x 10~ 2 , 10~ 2 , 10~ 3 , 10~ 4 , 10~ 5 and 10~ 6 M ) and all 
other parameters remaining fixed. It is important to stress that 
this series of models is not intended to reflect an evolutionary 
sequence as we do not change the dust properties accordingly 
(e.g. dust grain sizes, dust-to-gas mass ratio, settling). The goal 
is to explore a range of physical conditions, study their impact 
on the disk chemistry and analyze how this impacts the cooling 
radiation that will be probed e.g. with the Herschel satellite. 
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Fig. 2. From top to bottom row: disk models with 1CT 2 , 10~ 3 , 1CT 4 and 10~ 5 M Q . The first column shows the total hydrogen number 
density logn^ H > with white contours indicating gas temperatures of 100, 300, and 1000 K. The second column shows the UV 
radiation field strength^ (91.2 - 205.0 nm) from full 2D radiative transfer. The third column shows the gas temperature log T g with 
white contours indicating dust temperatures of 20, 40, 100, 300 and 1000 K. 
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Fig. 3. From top to bottom row: disk models with 10 2 , 1CT 3 , 10 4 and 1CT 5 M . The first column shows the C + density with white 
contours indicating the critical density for the [Cn] line and densities of 10 6 and 10 8 cirT 3 and a red contour line where the PDR 
parameter x/ n (H) = 0.01. The second column shows the O density with white contours indicating the critical densities for the two 
oxygen lines, «<h) = 6 10 4 and 5 10 5 cm -3 , and a density of 10 8 cirT 3 . Extinctions of Ay = 0.1 and 1 are denoted by the red contour 
lines. The third column shows the CO abundance with white contours indicating gas temperatures of 50, 100, and 200 K and the 
blue contours indicating dust temperatures of 50, 100, and 200 K). 
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Fig. 4. Departure coefficient logn NLTE /n LTE for Oi in the 10 2 M disk model. Blue contours denote the total hydrogen number 
densities with logarithmic values annotated in units of cirT 3 . 





Fig. 5. Departure coefficient log hnlte/^lte f° r Oi in the 10 4 M disk model. Blue contours denote the total hydrogen number 
densities with logarithmic values annotated in units of cirT 3 . 



4. 1 . Physical structure 



Fig-Elprovides an overview of the computed physical quantities 
such as gas density, UV radiation field strength, optical extinc- 
tion, gas temperature and dust temperature resulting from our se- 
ries of Herbig disk models; each row presents a model with dif- 
ferent mass and shows three quantities, total hydrogen number 
density, UV radiation field strength, and gas temperature. The 
white dashed lines in the second column of this figure indicate 
theAy = 1 andAy = 10 surfaces. The 2.2xl0~ 2 M model is op- 
tically thick in the vertical direction out to ~ 400 AU. Reducing 
the disk mass by a factor 200 leaves only a dense ~ 10 AU wide 
optically thick ring and the low mass disk models (< 10~ 4 M ) 
have very low vertical extinction, Ay < 0.1, throughout the en- 
tire disk. 



The high mass models show two vertically puffed up re- 
gions, one around 0.5 AU (the classical "inner rim"), the other 
one more radially extended between 5 and 10 AU. The very low 
mass models (10~ 5 M and below) do not show the typical flar- 
ing disk structure as observed in optically thick disk models, but 
they are still very extended in the vertical direction. The extreme 
gas temperatures in these low mass models are a result of direct 
heating by the photoelectric effect on small dust grains (PAHs) 
and pumping of Fe n by the stellar radiation field. In these op- 
tically thin disks, coupling between gas and dust grains is not 
efficient. 



4.2. Chemical structure 

The particle densities of ionized carbon, oxygen and the abun- 
dance of CO are shown in Fig. [3] The oxygen abundance is 
rather constant throughout the disks (C/O ratio ~ 0.45), being 
only lower by a factor of two where the CO and OH abundances 
are high; OH forms at high abundance inside 30 AU above an 
Ay of a few, thus co-spatial with the warm (> 200 K) CO gas. 
Extremely low abundances occur only if significant amounts of 
water form close to the midplane inside 1-10 AU (2.2 10~ 2 - 
10~ 3 M models). 

Carbon is fully ionized in the disk models with masses be- 
low 10~ 4 M Q . The C + abundance is complementary to the CO 
abundance as the neutral carbon layer between them is very thin. 
The C + /C/CO transition is governed by PDR physics and can be 
well defined using the classical PDR parameter x/ n {n)- Above 
X/ n (n) = 0.01 (see Fig. 01, carbon is fully ionized and its den- 
sity structure resembles that of the total gas density (two puffed 
up inner rims). The column density of the ionized carbon layer 
is always smaller than a few times 10 17 cm -3 . The mass of C + 
in the irradiated layers of the disk is roughly constant until the 
optically thin disk limit is reached in which case it is given by 
M gas e(carbon) (mc/jumn) (see Table|4j>. 

High abundances of CO can be found down to disk masses of 
10~ 4 M . Below that, the entire disk becomes optically thin (ra- 
dially and vertically), thus reducing the CO abundances even in 
the midplane to values below 10~ 6 . In the models with masses 
above 10 3 M , densities in the inner most region are hig h 
enough to form a large water reservoir (Woitk eet al.L 12009b). 
At densities in excess of ~ 10 11 cm 4 , the water formation con- 
sumes all oxygen, thus limiting the amount of CO that can form 
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Fig. 6. Oi 3 Pi level population numbers from Amc, log«AMC, 
and escape probability (ES), lognpRo, for the 1CT 2 M disk 
model (a, b) and the 10~ 4 M disk model (c, d). Blue contours 
denote the total hydrogen number densities with logarithmic val- 
ues annotated in units of cm -3 . 



in the gas phase. In the upper disk layers, the CO abundance 
is very low due to the combined impact of UV irradiation from 
the inside (central star) and outside (diffuse interstellar radia- 
tion field). Fig. [3] also shows that the temperature of CO in the 
outer disk decouples from the dust temperature (white contours: 
gas, blue contours: dust), even though differences are generally 
within a factor two (r gas > 7d ust ). This is relevant for the CO low 
rotational lines that are predominantly formed in those regions 
(see Sect. E3- 

Table 4. Characteristics of selected species in the Herbig disk 
models. (T g ) and (7a) are the mass averaged gas and dust 
temperature' 1 '. 



M disk [Ms] 


Species 


Mass [M Q ] 


(T g ) [K] 


<7d> [K 




O 


4.2(-5) 


61 


34 


2.2 1(T 2 


c + 


1.3(-7) 


96 


46 




CO 


6.0(-5) 


62 


35 




o 


2.1(-5) 


61 


35 


io- 2 


c + 


1.1 (-7) 


99 


46 




CO 


2.7(-5) 


61 


35 




o 


2.5(-6) 


64 


34 


1CT 3 


c + 


9.7(-8) 


79 


34 




CO 


1.8(-6) 


79 


42 




o 


3.2(-7) 


77 


39 


irr 4 


c + 


4.9(-8) 


75 


37 




CO 


9.2(-8) 


103 


51 




o 


7.0(-8) 


92 


50 


io- 5 


c + 


2.0(-8) 


91 


48 




CO 


3.6(-12) 


576 


126 




o 


5.1(-8) 


53 


43 


itr 6 


c + 


1.7(-8) 


53 


43 




CO 


6.1(-15) 


89 


54 



(1): <T> = f ( '„~ ) "r''-)rfv V ' s P ec ' es density n sp and volume element dV. 



4.3. Self-similarity 

The six disk models show a large degree of self-similarity in 
their temperature and chemical structure. Lowering the disk 
mass removes mainly the thick midplane and the innermost 
dense regions. In that sense, this series of decreasing mass zooms 
in into disk layers further away from the star and at larger 
heights. The reason for this self-similarity lies partly in the two- 
direction escape probability used to compute the gas temperature 
and partly in the chemistry being independent of neighboring 
grid points (no diffusion or mixing). The solution in the disk is 
mostly described by the local xl n {n)- 

4.4. Strength of UV field 

To assess the impact of the UV field on the disk structure 
and line ratios, we computed two additional disk models with 
10~ 2 M disk masses, using effective temperatures of 9500 and 
10500 K. This range reflects the typical temperatures encoun- 
tered for Herbig Ae stars. To isolate the effect of UV irradiation, 
we keep the luminosity constant, so that a change in effective 
temperature changes only the fraction of UV versus optical irra- 
diation. 

Increasing the stellar effective temperature to 10500 K leads 
to a vertically more extended disk structure, thus pushing the 
H/H2 transition slightly outwards. Though the dust temperatures 
are unaffected (they depend rather on total luminosity), the mass 
averaged gas temperatures increase by up to a factor two for cer- 
tain species such as CO and O. The change for C + being in the 
uppermost tenuous surface is more dramatic; its mass averaged 
temperatures increase from ~ 100 K (T e s = 8460 K, Tabled to 
330K(r eff = 10500 K). 

4.5. Dust opacities 

In a similar way, we varied the dust opacity in the 10~ 2 M© 
disk mass model, using first very small grains a m ; n = 0.05 
to a max = 1 micron and then only large grains a m i n = 1 to 
flmax = 200 micron. They represent the two extremes of grain 
size distribution ranging from rather pristine ISM grains to con- 
ditions appropriate for more evolved dust in very old disks. 

Dust opacities impact disk physics in two ways. First, in- 
creasing the average grain size decreases the opacity and thus 
the optical depth in the models. The dust temperature de- 
creases mostly, except for the midplane regions inside 100 AU. 
However, the second — more important — effect is a decrease 
of effective grain surface area, thus decreasing the efficiency of 
gas-dust collisional coupling. As a consequence, the gas temper- 
ature decouples from the large grains even in the disk midplane, 
leading to higher gas temperatures everywhere and thus a verti- 
cally more extended disk. Given the fact that gas dust coupling 
is no longer efficient, the initial assumption that gas and dust are 
homogeneously mixed would have to be revisited. If most heat- 
ing is provided by PAHs, this is not an issue as those will stay 
well mixed with the gas. 

5. Line emission 

In the following, we performed radiative transfer calculations 
on the grid of Herbig Ae disk models to understand the spatial 
and physical origin of the two Gasps tracers [Cn], [Oi], and the 
frequently observed sub-mm lines of CO. Besides a wealth of 
published CO observations of protoplanetary disks to compare 
to and test the physics and chemical networks of our models, the 
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Table 5. Results for the fine structure emission lines from various disk models and test runs. In the model column, MC, LTE and 
ES denote the use of level populations from Monte Carlo line radiative transfer, LTE and escape probability, respectively. Line and 
continuum fluxes are shown for a distance of 131 pc and an inclination of 45 degrees. Numbers in parentheses indicate powers of 
ten, i.e. 1.93(3) denotes 1.93 x 10 3 . The peak separation is defined as the distance between the two maxima in the double peaked 
line profile measured in km/s. 









O i 63 yUm 






M t 1 nm 
\ J 1 L<iD yum 






II I J / . / fJlTl 








1 line 


F . 

1 cont 


P^qI^ cpn 
L CtlJv SCU. 


FY 
1 line 


F 

' cont 


r CtiJv acu. 


FY 
1 line 


F 

1 cont 


Pf*il^ cpn 

L CtllV 3CU. 






no -18 W/m 2 l 






no -18 W/m 2 

|_±VJ VV/lll 




rtm/<;1 
Liviii/ & j 


no~ 18 W/m 






2.2 1(T 2 


MC 


1.93(3) 


1.44(2) 


1.9 


2.13(2) 


1.25(2) 


2.6 


5.39(1) 


1.14(2) 


2.2 


10 2 


MC 


1.29(3) 


1.03(2) 


2.4 


1.18(2) 


6.57(1) 


3.0 


4.69(1) 


5.91(1) 


2.2 


10 3 


MC 


6.57(2) 


2.58(1) 


2.5 


3.89(1) 


8.84 


5.9 


2.98(1) 


7.69 


2.2 


10 -4 


MC 


2.65(2) 


4.56 


5.9 


8.06 


9.09(-l) 


6.2 


1.40(1) 


1.01 


2.4 


10 


MC 


4.19(1) 


9.89(-l) 


7.7 


1.11 


1.72(-1) 


6.0 


4.53 


2.71(-1) 


2.3 


10~ 6 


MC 


1.99 


1.50(-1) 


12.2 


5 (-2) 


3(-2) 


10.3 


2.77 


L43(-l) 


2.2 


i n-'l 
10 


LTE 


2.25(3) 


1.05(2) 


2.2 


1.30(2) 


6.63(1) 


2.7 


5.22(1) 


5.96(1) 


2.2 


10 


ES 


1.27(3) 


1.05(2) 


2.3 


1.23(2) 


6.63(1) 


2.9 


4.93(1) 


5.96(1) 


2.2 


io~ 2 


no O i > 2000 K 


1.29(3) 


1.03(2) 


2.4 


1.18(2) 


6.57(1) 


3.0 








io- 2 


«i„. ou t = 0.5, 500 AU 


1.23(3) 


1.06(2) 


2.6 


1.14(2) 


6.16(1) 


3.2 


2.26(1) 


5.41(1) 


2.7 


io- 2 


fiin.01,1 = 0.5, 300 AU 


9.40(2) 


1.01(2) 


3.3 


8.99(1) 


4.98(1) 


4.0 


7.33 


4.66(1) 


3.5 


10- 2 


fiin.o Ut = 0.5, 100 AU 


5.74(2) 


8.27(1) 


5.6 


4.18(1) 


3.19(1) 


6.0 


1.35 


2.76(1) 


6.3 


io- 2 


«in.ou. = 0.5,30AU 


2.40(2) 


5.03(1) 


11.2 


1.26(1) 


1.47(1) 


11.4 








10- 2 


flin,out = 30, 100 AU 


3.70(2) 


3.42(1) 


5.5 


3.08(1) 


1.75(1) 


6.1 








io- 2 


«in.out = 30,700AU 


1.11(3) 


6.13(1) 


2.5 


1.08(2) 


5.29(1) 


3.0 


4.64(1) 


4.82(1) 


2.2 


10- 2 


Rm,om = 100, 700 AU 














4.59(1) 


3.31(1) 


2.2 


IO- 2 


r eff = 9500 K 


5.53(3) 


1.68(2) 


2.8 


4.56(2) 


8.60(1) 


2.9 


1.50(2) 


7.63(1) 


2.3 


io- 2 


T e g = 10500 K 


9.03(3) 


1.78(2) 


3.1 


6.75(2) 


8.85(1) 


2.8 


3.17(2) 


7.84(1) 


2.3 


IO' 2 


amm.max = 0.01, 1 /HID 


6.48(2) 


1.42(2) 


2.3 


3.24(1) 


3.53(1) 


2.5 


6.30(1) 


2.85(1) 


2.2 


io- 2 


a m m.max = 1,200/im 


3.22(3) 


1.44(2) 


2.7 


3.17(2) 


7.21(1) 


3.2 


2.39(1) 


6.42(1) 


2.2 



low rotational CO lines are part of ancillary projects to comple- 
ment the Herschel Gasps project with tracers of the outer cold 
gas component in disks. 



5.1. [Cll]158fim line 

The fine structure line of ionized carbon arises in the outer sur- 
face layer of the disk. For values of x/n(H) smaller than 0.01, 
carbon turns atomic/molecular (see Fig. [3] left column). This 
limits the C + column density and hence the [Cn] 158 //m line 
emission. Except in the very inner disk, r < 1 AU, the line never 
becomes optically thick. 



5.1 .1 . Line formation regions 

The line can be easily excited (E u = 91.2 K) even in the disk 
surface at the outer edge; hence the total [C n] emission from the 
disk is dominated by the 100-700 AU range and probes the gas 
temperature in those regions (r gas + 7d us t)- At these distances, 
the column density of C + decreases with disk mass leading to 
a potential correlation of the [Cn] 158 /mi line emission with 
total disk gas mass. However, as shown in Table [5] the total line 
emission depends sensitively on the outer disk radius. The total 
[C n] line emission is thus degenerate for disk mass and outer 
radius. 

In addition, the surrounding remnant molecular cloud mate- 
rial will also emit in the [C n] line, the only difference being in 
general lower densities and temperatures than those encountered 
in our protoplanetary disk models. The mass averaged gas tem- 
perature of C + in the disk is ~ 90 K. Densities range from up to 
10 5 cirT 3 in the outer disk (700 AU) to several times 10 s cm -3 
in the regions inside 10 AU (close to the;^/«<H>=0.01 layer). 



5.1.2. LTE versus escape probability versus Monte Carlo 

Due to the low critical density of this line, the emission forms 
largely under LTE conditions. Deviations from LTE are small, 
less than 10%, and grow towards lower disk mass models (10~ 4 
- 20%, 10~ 5 M Q - 40%). Escape probability and Monte Carlo 
line fluxes agree well within 5-10%, the typical uncertainty that 
can be expected from the re-gridding. Also here, the larger dis- 
crepancies are found in the lower mass disk models. 

5.2. [01]63 and 145^m lines 

Since our disk models span a much wider range of temperatures 
and densities than those found in molecular clouds and shocks, 
we encounter in this paper also different regimes for the forma- 
tion of these two fine structure lines. The following paragraphs 
explore this in more detail. 

5.2.1. Line formation regions 

The column densities at which the 63 fim line becomes optically 
thick can be approximated («/ ~ no, n u ~ 0) as 

_ g, 8^v 3 Av D 

-Wthick T— (11) 

gu A ld C^ 

where g„ and gi are the statistical weights of the upper and 
lower level, A u \ the Einstein A coefficient for the respective 
line transition with the frequency v. As an approximation, the 
Doppler width Avd is assumed to be 1 km/s (within ProDiMo, 
the Doppler width is calculated from the actual sum of thermal 
and turbulent broadening). The same formula can be used for the 
145 jum line, if we take a factor into account correcting for the 
true level populations. 

The surface of the inner 30 AU of the 10~ 2 M disk model 
is very hot with gas temperatures of several thousand K and typ- 
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ical densities above 10 7 cm 3 . Thus, the relative level popula- 
tions with respect to the total oxygen number density no, in LTE, 
follow from the ratios of the statistical weights, «2/«o = 0.56, 
«i/«o = 0.33, «o/«o = 0.11 (see Table [3] for the notation — 
nj with g = 27+1 being the statistical weight of the level). 
Under these physical circumstances, both O i lines become opti- 
cally thick at column densities of ~ 3 10 17 cm 2 . Such column 
densities are reached inside 30 AU, even in our lowest mass disk 
model. However, the contribution to the integrated line emission 
from this hot gas is negligible. Removing the contribution from 
hot gas by setting the O i abundance to zero for gas temperatures 
in excess of 2000 K, results in line flux changes that are less than 
1 %. 

In the 30-100 AU range, gas temperatures are much lower 
(few hundred K) and only about 10 % (3 %) of the oxygen atoms 
reside in the upper level of the 63 (145) /mi line. For those re- 
gions, the 145 //m line becomes optically thin in the 10~ 3 and 
10~ 4 M disk models. Emission from this ring dominates the to- 
tal line emission. We will get back to this point later. Models 
that include X-ray hea ting and ionization (Nom ura et al.L 2007; 
Meii erink et all 120081) indicate that the impact of those addi- 
tional processes is mostly relevant inside 30 AU. Beyond that, 
UV processes dominate the gas energy balance and chemistry. 
Hence, we do not expect X-rays to significantly alter the [O i] 
line emission for disks with masses > 10~ 4 M . In lower mass 
disks, the [Oi] emission originates closer to the star and here 
X-rays could have an impact by increasing the gas temperatures. 

The outer disk, beyond 100 AU, is fairly cold (T g < 200 K). 
Here, level population numbers of the J=0 and J=l level are ex- 
tremely small (below 0.5 %). The peak separation 6v — 1 .9 km/s 
of the [Oi]63 /mi line in the 2.2 x 10~ 2 M disk model (see 
Fig. 13 indicates that the emission comes mostly from inside 
380 AU; the [Oi] 145 /mi originates slightly closer to the star, 
inside 200 AU (<5v = 2.6 km/s). The disk surface outside of 
100 AU accounts for roughly 50 % of the O i emission. For the 
lowest mass disk model, the physical conditions in this region 
such as gas/dust temperature and densities get close to molecu- 
lar cloud values of 7g = 50 K, T d = 40 K, n <H > = 10 5 cm -3 . 

In a series of radiative transfer calculations for the chemo- 
physical structure from the 10~ 2 M disk model, we studied 
the fraction of the emission coming from regions with densi- 
ties above «<h>- Varying that density from the critical density 
of the 63 /mi line, «<h> = «crit = 5 10 5 cirT 3 , to a density of 
10 8 cm 3 , shows that the regions below the critical density do 
not contribute very much to the total line emission. The emis- 
sion gradually builds up with increasing density, showing that 
the line originates from regions up to an extinction of Ay ~ 0.1. 

5.2.2. LTE versus escape probability versus Monte Carlo 

LTE level population numbers systematically overestimate line 
fluxes by up to 70%. In the regions where most of the Oi emis- 
sion arises, the upper levels are significantly less populated than 
in LTE while the ground state is overpopulated with respect to 
LTE (Fig. |4] and [5). In the regions inside 30 AU, the H2 abun- 
dance is low, thus the main collision partners are atomic hydro- 
gen and electrons. Outside 30 AU, hydrogen is predominantly in 
molecular form. 

The level population numbers computed with the simple 
escape probability assumption in ProDiMo yield only slightly 
higher line fluxes, within 10 % for the 1 lines and ~ 5% higher 
continuum fluxes. Calculating the line emission from Monte- 
Carlo radiative transfer requires a re-gridding of the ProDiMo 



model results. Fig.|6]shows the level population numbers of the 
3 Pi level of oxygen (upper level of the 63 /mi line) for two mod- 
els. Differences arise mainly in the very optically thin regions 
either at the disk surfaces or near the outer edge. In these areas, 
the maximum escape probability following our two directional 
escape probability approach is 0.5, while it should be rather 1.0 
in very optically thin environments. Hence, in those areas the 
Monte Carlo approach gives a better result. 

Since continuum differences should largely be due to the re- 
gridding, the estimated error stemming from the interpolation 
onto a different grid is of the order of 5 %. Considering that, the 
intrinsic difference in line fluxes of both line radiative transfer 
methods is fairly small. The visible differences in the level pop- 
ulation numbers are largest in the regions that do not contribute 
significantly to the total line emission (well above the critical 
density iso-contour in Fig. |6), i.e. the very tenuous surface and 
outer (r » 100 AU) regions of the disk. 

5.3. CO sub-mm lines 

Line fluxes have been calculated for the first three molecular 
transitions in CO, across the full mass range of disk models. The 
line fluxes and continuum fluxes are plotted as a function of disk 
mass in Fig. [8] 

The J=l-0, J=2-l and J=3-2 lines all show similar behaviour 
with disk mass, with the line fluxes initially increasing sharply 
with mass before levelling off for Mdisk > 10~ 4 M . This similar- 
ity results in largely uniform line ratios, with a spike at 10~ 5 M Q . 
The continuum varies almost linearly with disk (and hence dust) 
mass, although with slightly more emission from the lower mass 
models than a strict linear relationship would give. This is due 
to the optically thin low mass disks having more penetrating UV 
dust heating in the midplane than the higher mass optically thick 
disks. 

The calulated line fluxes all exhibit a jump of four orders of 
magnitude on moving from the 10~ 5 to the 10 4 M model. This 
discontinuity is also seen in the line ratios, with all three exhibit- 
ing a significantly higher ratio for the 10~ 5 M model than for 
the others. The sudden drop in emission below 1O~ 4 M corre- 
sponds to a fall-off in Ay. There is suddenly almost no region 
in the disk with Ay > 0.1, resulting in a maximum CO abun- 
dance of less than 10~ 6 compared with ~ 10 -4 for the higher 
mass models. The small remaining region with Ay > 0.1 is quite 
hot, with gas temperatures around 500 K, giving a spike in the 
line ratios at 10~ 5 M . 

5.3.1. Line formation regions 

The three CO lines form at an intermediate height in the disk, 
between the warm upper layer and the cold midplane (see also 
Ivan Zadelhoff et aU |2001|) and are generally optically thick for 
disk masses down to 10~ 4 M . The submm line formation re- 
gion is radially very extended, with significant contributions to 
the total line flux from the entire disk. The results of an experi- 
ment varying the outer radius of the disk region sampled in the 
re-gridding procedure are plotted in Fig. [9] The total J=3-2 line 
flux is seen to vary linearly with R out indicating that the line orig- 
inates from the full radial extent of the disk. The same behaviour 
is seen in the J=l-0 and J=2-l lines. The linear trend is caused 
by a combination of the radial gas temperature gradient in the 
CO emitting layers and surface area. The continuum shows a 
slightly different behaviour, with a greater proportion of the in- 
tegrated flux from outer radii. 
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Fig. 7. Top row: [Oi] 63/mi line emission flux from the 1CT 2 M and 10~ 4 M model at a distance of 131 pc and an inclination of 
45 degrees. The inserted image shows the continuum subtracted integrated logarithmic line intensity as a function of sky position. 
The angular size of 16" corresponds to 2100 AU at the distance of 131 pc. Bottom row: [Cn] 158/im line emission flux from the 
10~ 2 M model (note the different velocity range) and [Oi] 63 //m line emission flux from the 10~ 5 M Q model. The color scale of 
the inserted image is the same in all panels (/ max = 3 10~ n , 7 m j n = 3 10~ 14 W/m 2 /Hz/sr). 



The line profiles for the three transitions are generally very 
similar for a given disk model, and indeed across the computed 
mass range. Narrow peak separations (6\ =1-2 km/s) indicate 
that the emission is coming from the entire disk inside ~ 700 AU 
and will thus be dominated by the outer regions that contain 
more surface area. 



5.3.2. LTE versus escape probability versus Monte Carlo 

The continuum results obtained without re-gridding (ES) are 
in good agreement with those from the Monte Carlo radiative 
transfer, within 3% in all cases. This indicates that re-gridding 
does not present an issue when considering sub-mm fluxes. The 
escape probability line fluxes deviate slightly from the Monte 
Carlo and LTE fluxes for the low disk masses, up to ~ 50% lower. 



This indicates again the limits of the two direction escape proba- 
bility approach for the optically thin case of very tenuous disks. 



The escape probability fluxes are however in good agreement 
(within 3%) for disk models with masses larger fhanl0~ 4 M o . 
The LTE line fluxes are also in good agreement (within ~ 1 %) 
for the entire mass range (see Fig. [Tol l, indicating that LTE is 
a valid approximation for these low CO lines. However, r eas = 
Tdust is not a valid approximation for this set of disk models, be- 
cause the gas temperatures in the outer regions, where the CO 
lines arise, deviate by up to a factor two from the dust tempera- 
tures (see Table |4]for mass averaged gas and dust temperatures 
of CO). 
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Fig. 8. The top row shows integrated line fluxes as a function of disk mass, for (from left to right) the CO J=l-0, J=2-l and J=3-2 
rotational lines. The second row shows the continuum fluxes at the corresponding line center wavelengths as a function of disk mass. 
The third row shows the CO line ratios as a function of disk mass:(from left to right) J=2-l/l-0, J=3-2/l-0, J=3-2/2-l. Dashed lines 
indicate computed line ratios for optically thick lines in LTE. 
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5.3.3. CO line ratios 

To understand the values of the CO line ratios for high mass 
disks in Fig. [8] we start with a number of assumptions: 1) The 
continuum is small 7 cont <K / v ; 2) The CO line forms under opti- 
cally thick LTE conditions with a universal line profile 0(v) such 
that I v = B v (7gas) 0( y ) where B v is the Planck function; 3) The 
temperature r gas is constant throughout the line forming region. 
The line flux can then be expressed as 

^line = jj (/y - /cent) dCl dv * A B v (T gas ) J 0(v) dv , (12) 

where A is the disk surface area as seen by the observer. 

If we assume a square line profile function in velocity space 
such that cp(v) — 1 if -4r < v < +4? and cf>(v) = otherwise, we 
can re-write equation ([121 in the Rayleigh-Jeans limit as 



v3 



(13) 



Fig. 9. Integrated CO J=3-2 line flux for the 10 J M Q model, plot- F^ e =AB v (T gss ) Av=A2kT m Av- 
ted against the outer radius sampled by re-gridding. c 

using — = For line ratios, the quantities A, r gas and Av are 
identical, and hence 



■Fline.l 
fline,2 



(14) 
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Fig. 10. Relative difference in percent between Monte Carlo CO 
line fluxes and CO line fluxes calculated using the escape prob- 
abilty assumption (solid line) and between Monte Carlo CO line 
fluxes and LTE CO line fluxes (dotted line); results are shown as 
a function of disk mass 



The corresponding values for the three CO line ratios are over- 
plotted with dashed lines on Fig. [8] and agree indeed well with 
the limiting behaviour of the line ratios at high disk masses, 
A^disk ^ 1CT 3 . Thus, we conclude that CO line formation in high 
mass disks occurs under optically thick LTE conditions. At lower 
masses the line ratios deviate from these values since the as- 
sumption of optical thick LTE no longer holds. 

Changing the effective temperature of the central star or 
changing the dust properties affects the CO lines to a lesser ex- 
tent than the fine structure lines. Changes to the continuum and 
the lines are generally within 50%. Even though the formation 
height of CO might change slightly, the excitation conditions are 
still close to LTE. While the CO mass in the models stays within 
a few percent of the original model (except in the ISM grain 
model, where it is 26% smaller) the CO mass averaged temper- 
ature changes up to a factor two. The effect on the line fluxes 
correlates directly with the temperature change, with the higher 
temperatures giving systematically higher fluxes. At these radio 
wavelengths, the change in continuum emission amounts to less 
than 25%, even between the models with different dust proper- 
ties. The most extreme model is the one with an effective tem- 
perature of 10500 K where the mass averaged CO temperature 
rises to 1 10 K and hence the line fluxes increase by ~ 60 - 70%. 

6. Discussion 

The goal of this study is a first assessment of the diagnostic 
power of the fine structure lines of neutral oxygen, ionized car- 
bon and of the CO submm lines in protoplanetary disks. This 
work prepares the ground for a systematic analysis of a much 
larger and more complete grid of protoplanetary disk models. 
Based on the detailed radiative transfer calculations of the pre- 
vious section, we discuss here the potential of the various lines 
in the context of deriving physical disk properties such as gas 
mass, gas temperatures and disk extension. 

6.1. [Oi] 63 line 

Using o ne dimensional mod eling of photon-dominated regions 
(PDRs), iLiseau et all d2006h find that the total intensity of the 
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Fig. 11. [Oi] 63/145 fim line ratio as a function of disk gas mass. 
In the high disk mass regime, both lines are optically thick and 
the line ratio levels off around 10. At very low disk masses both 
lines are optically thin and the line ratio of ~ 40 corresponds 
to ga s temperatures above 1 00 K and a wide range of densities 
(see iTielens & Hollenbachl |l985). Plus signs denote the stan- 
dard mass sequence, open circles models with varying outer ra- 
dius (100, 300, 500 and 700 AU). The solid box show the range 
of line ratios that occurs with varying dust properties to the ex- 
tremes, ISM type grains (a m i n = 0.05, a max = 1 fim) and large 
grains (a min = l,a max = 200 fan). 
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Fig. 12. [Oi] 63/[Cii] 158 fim line ratio as a function of disk gas 
mass. Plus signs denote the standard mass sequence, open circles 
models with varying outer radius (500 and 700 AU). 



63 fim line can be used as an indicator of the gas temperature in 
outflows from young stellar objects. From detailed two dimen- 
sional modeling of protoplanetary disks, Woit ke et al.l (|2009a) 
conclude that the [O i] 63 //m line alone will indeed provide a 
useful tool to deduce the gas temperature in the surface lay- 
ers of protoplanetary disks, especially the hot inner disk sur- 
face. However, spatially unresolved observations that measure 
only integrated [O i] 63 fim line fluxes in sources at typical dis- 
tances of 140 pc will not be able to detect the presence of hot 
gas very close to the star (inside ~ 10 AU). The total line flux 
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Fig. 13. [Oi] 63/145 fjm versus [Oi]63/[Cn] 158 /mi line ratio 
in the series of Herbig Ae disk models. Plus signs denote the 
standard mass sequence, open circles models with varying outer 
radius (500 and 700 AU). The data points for smaller outer radii 
fall well outside the plotting range as indicated by the arrow (see 
Table[5]l. The star denotes the location of the 10~ 2 M disk model 
with small grains (a m j n = 0.05, a max = 1 ^m) 



is dominated by emission coming from the cooler outer regions 
(TableEJ. 



6.2. [Oi] 63/145 line ratio 

The diagnostic power of the [Oil 63/145 /m i line ratio has been 
discussed by Tiel ens & Hollenbachl d 1985b for PDRs and by 



by llieler 

iLiseau et al.l ([2006) for outflows around young stellar objects 
(YSOs). Very low ratios of a f ew can be explained if both li nes 
are optically thick (see Fig. 2 o flTielens & Hollenbachl fl985h . In 
that case the diagn ostic power of this line ratio is very limited. 
ILiseau et al.l ((2006) find that even in the optically thin limit, the 
line ratio can be degenerate with low temperature molecular gas 
giving the same result as high temperature atomic gas. This is 
due to the different colli sional cross sections fo r O-H and O-H2. 
From ISO observations. lLorenzetti et all d2002l) found ratios be- 
tween 2 and 10 for Herbig Ae/Be stars. They derive from PDR 
modeling of a clumpy medium that the low line ratios of a few 
would require very high column densities of N(n) ~ 10 22 cm -3 
(Ay ~ 18-27 mag) and conclude that better models of the emit- 
ting region are required for the physical interpretation of this 
data. Studyin g the outflows of several hundred YSOs with ISO, 



lying the o 

ILiseau et all d2006l) found [O I] 63/145 /mi line ratios between 
~ 1 and 50. They suggest that confusion with relatively cool but 
optically thick envelope gas (T < 100 K) could explain the low 
line ratio in some cases. 

Fig. fTTI shows that the [Oi] 63/145 line ratio increases to- 
wards lower masses (10 for ~ 10~ 2 Mq to 40 for M disk < 
10~ 4 M ). Our line ratio for massive disks corresponds well to 
the expected value for both lines being optically thick and aris- 
ing from relatively warm gas, T > 100 K, that resides in the 
surface layers of the disk. As the disk mass decreases, first the 
145 and then the 63 fim line become optically thin. Eventually, 
the optically thin limit is reached for disk masses smaller than 
10~ 4 M , leading to line ratios of ~ 40. In a simple one dimen- 
sional model, this corresponds to gas temperatures in excess of 



100 K and a wide range of densities. The line profiles show an 
increasing peak separation towards lower mass disks (see Fig.[7]l. 
In fact, the line emission in low mass disks is dominated by disk 
material close to the star (r < 100 AU), while in high mass disks, 
the emission comes from the disk surface out to ~ 200 AU and 
is hence dominated by the outer regions (30 < r < 200 AU). The 
sequence of the 10~ 2 M disk model with varying outer radii 
(dotted line in Fig. fTTb shows that the [Oi] 63/145 line ratio is 
indeed not sensitive to the disk outer radius unless the disk starts 
to become smaller than the region where most emission origi- 
nates (see the R out = 100 AU model). 

The effective temperature and thus amount of UV radiation 
from the central star has a negligible impact on the line ratio in 
the temperature range of Herbig Ae stars. Both lines increase by 
the same amount (up to a factor 7 for the individual line fluxes 
for T e ff = 10500 K compared to the standard case of T e s = 
8460 K) due to the overall warmer gas temperatures (the oxygen 
mass does not change by more than 5%). The change in line ratio 
due to a factor 10 decrease in disk mass is much larger than the 
change due to a different effective temperature of the star. 

On the other hand, the grain size distribution does have a 
larger effect on the line ratios as expected from the impact on 
the disk structure (see Sect. |4.5I >. The line ratio is 10% smaller 
for the disk model with large micron sized dust grains. Using 
only small sub-micron sized grains (similar to the ISM compo- 
sition) increases the line ratio by a factor two, thus mimicking a 
disk mass that is lower by more than a factor ten. Since the disk 
with small ISM grains is flatter and reaches thus higher particle 
densities, O depletes more readily into species such as water and 
water ice, hence decreasing the overall vertical column density 
of neutral oxygen (making the lines thus marginally optically 
thin). However, the grain size distributions chosen here present 
the extreme ends and are easily distinguished from looking at 
theSED. 

The typical hydrogen number densities in the emitting region 
range up to 10 s crrT 3 , higher than those generally covered in 
PDR models (even in clumpy PDRs). The vertical column den- 
sity of neutral oxygen reaches values of 10 20 cm~ 3 and tempera- 
tures above 200 K inside 10 AU and 10 19 cirT 3 and temperatures 
below 100 K in the outer disk (10~ 2 M disk model). The large 
difference to PDR con ditions can expla i n the p roblems and high 
extinction values that lLorenzetti et al.l (120021) found from their 
classical PDR analysis. However, our highest mass disk mod- 
els do not show line rat ios as low as 2, the most extreme cases 
of the YSOs studied by ILiseau et al. (2006). For embedded ob- 
jects, confusion with relatively cool but optically thick envelope 
gas (T < 100 K) is still a viable explanation of the low ratios. 



6.3. [Cn]158ijm line 

The [Cn] 158 /mi line emission is a very sensitive tracer of the 
outer disk radius. We come back to this point in Sect. 16.51 As 



expected, the integrated line flux from our models scales as 7? 2 ut . 
This makes it more difficult to detect smaller size disks that are 
only of the order of 100 AU in size. 

Due to its low excitation temperature, this line is also 
strongly affected by confusion with any type of diffuse and 
molecular cloud material in the line of sight. The problem can 
be addressed with Herschel/PACS by using the surrounding pix- 
els on the detector (5x5 pixels with the disk being spatially 
unresolved) for a reliable off-source flux determination. 

Changing the dust properties to very large grains results in a 
much lower [Cn] 158 fim line flux, the reason being that above 
the x/ n (H) = 0.01 layer outside of 100 AU, the gas tempera- 
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ture drops below the excitation temperature for this line. The 
gas inside 30 AU becomes much warmer due to the gas-dust de- 
coupling, thereby generating a more pronounced secondary rim. 
This casts an efficient shadow on the outer disk, thereby enhanc- 
ing the CO formation (H2 and CO self-shielding) and thus the 
molecular line cooling. 

6.4. [Oi]63/[Cn] 158 line ratio 

The [0 1] 63/[C 11] 158 ratio depends even more on disk mass than 
the [Oi] 63/145 line ratio. It would thereby be a prime diagnos- 
tic for relative disk mass estimates. However, while the ratio 
of the two oxygen lines does not change very much with disk 
outer radius, the [0 1] 63/[C 11] 158 ratio does depend strongly on 
it (Fig. ITSb . The 10~ 2 M Q disk model with an outer radius of 
100 AU has a more than 10 times higher [Oi] 63/[Cn] 158 line 
ratio than the same model with an outer radius of 700 AU. 

Because of this degeneracy with disk size, the diagnostic 
value of this ratio alone is limited for individual objects. Its use 
would require the a-priori knowledge of many other disk param- 
eters such as the disk outer radius, grain sizes and composition. 
Rather we suggest that it could be used in a statistical way on 
a large sample of disks that share the same classification group, 
e.g. based on their IR SEDs. Alternatively, a combination of the 
[Oi] 63/145 and [Oi]63/[Cn] 158 line ratios could be used to 
break the degeneracy with disk size. Fig. [T3l shows clearly that 
disks with a smaller outer radius, e.g. 300 AU, lie well to the 
right of the standard relation, while the presence of an inner hole 
for example (30 AU size here) does not affect the location in this 
diagnostic diagram. 

As for the oxygen fine structure line ratio, the 
[0 1] 63/[C 11] 158 ratio depends strongly on the chosen grain size 
distribution. Additional SED observations are required to con- 
strain the dust parameters to a reasonable range. It is important 
to note that the change in line ratios in the [Oi] 63/145 versus 
[Oi]63/[Cii] 158 plot has a direction different from what we 
expect from a change in disk mass. Also, the [Oi] 63/[Cn] 158 
ratio is less sensitive to the effective temperature of the star than 
it is to disk mass. 



6.5. CO submm lines and line ratios 

The low excitation CO rotational lines are often used to prove the 
existence of primordial gas in protoplanetary disk and in combi- 
nation with opticall y thin line of isotopologues to me asure the 
total gas mass (e.g. IThi et all 1200 it IPanic et all 12008b . Due to 
the high Einstein A coefficients and high abundances of 12 CO 
in molecular regions, the low rotational CO lines become opti- 
cally thick at low Ay. Thus, they originate m ainly in the optically 
thin disk surface (Ivan Zadelhoff et all 1200 ll) . where gas and dust 
temperature decouple according to our models. However, as in 
the case of the [C n] line, CO submm lines can be contaminated 
by low density cold remnant cloud material. The amount of con- 
tamination can be estimated either from offset positions for sin- 
gle dish observations or through interferometry. 

Isotopologues of CO ha ve much lower abun dances of 
12 CO/ 13 CO = 77 in the ISM dWilson & Roodl [l994h and even 
higher values of 10 in disks (selective photodissociation; see for 
most recent work ISmith etaO. l2009t IWoods & Willacvi 120091; 
IVisser et"ail [2009). Thus, they can probe deeper into the disk 
and are — in the absen ce of freeze - out — excellent tools for 
measuring the disk mass. iPanicet al.l(l2008l) use for example the 



optically thin I3 CO J=2-l line from the disk around the Herbig 
star HD 169142 to infer a disk mass of 0.6 - 3.0 10~ 2 M . 

Our analysis shows indeed a relatively flat relationship for 
the 12 CO lines with disk gas mass, confirming thus that the low 
rotational CO lines by themselves are not a good tracer of gas 
mass (Fig. |8]l. The lines become optically thick even at very low 
disk masses and hence a reliable inversion of the line fluxes to in- 
fer gas mass proves impossible. Similarly there are degeneracies 
in the computed line ratios J=2- 1/1-0, J=3-2/l-0 and J=3-2/2-l 
which suggest limitations in the diagnostic power of these line 
ratios. 

The CO lines can be used to measure independently the outer 
disk radius (Sect. 15.3. It and thus mitigate th e problem of em- 
ployi ng the fC nl 158 line f or disk diagnostics. iDent et alj (120051) 
and P anic & HogerhejjdeJ {2009) demonstrate the power of sin- 
gle dish CO observations in deriving disk outer radii and incli- 
nations, thus constraining SED modeling that is rather insen- 
sitive to R out and deg enerate for inclination and inner radius. 
Hughes et al. (2008b|) note that the outer radius derived from 
dust continuum observations is always smaller than that de- 
rived from CO emission lines. Even though they can explain this 
within homogeneous gas and dust models (R out (gas) = /? out (dust) 
and a realistic outer edge, i.e. an exponential density profile), this 
finding reveals the problems in mixing quantities deduced from 
dust and gas observations. Since C + and CO arise both from the 
gas in the outer disk, we do not need to rely on dust observa- 
tions, but can instead use a gas tracer to measure the gas outer 
disk radius. 

The [C n] line probes in general slightly higher layers than 
the CO lines; hence the model results suggest that we can trace 
the vertical gas temperature gradient in disks. Table|4]shows that 
for the most massive disks, (7c+) is 50% larger than {Tco)- 

6.6. Comparison with observations for MWC480 

Since we used the Herbig star MWC480 to motivate our choice 
of parameters, we briefly discuss the line emission from the 
models together with observations available from the literature. 
The goal is not to obtain a best-fitting model, but to provide 
an impression of the model results within the observational 
context. Table [6] gives an overview of selected observational 
data for MWC480 (detections and upper limits) from the liter- 
ature. The CO mass derived from 13 CO 3-2 is 1.7 10~ 4 M and 
the total gas mass derived from the 1.3 mm continuum flux is 
2.2 ± 1.0 10" 2 M Q (IThi et all l200lh . The 12 CO and 13 CO lines 
sugge st an outer disk ra dius of ~ 750 and ~ 500 AU, respec- 
tively dPietuet all 120071) . 

The oxyg en fine structure lines are n ot detected with the 
ISO satellite dCreech-Eakman et all |2002|) . Assuming the ISO 
instrumental resolution as FWHM of a potential line (0.3 /mi for 
A < 90 /mi, 0.6 /mi for A > 90 /mi), the ISO 3<x upper limits for 
the [0 1] lines can be derived from 

Supper limit * 3<r V05^ FWHM W/m 2 . (15) 

This gives values of 4.14 10~ 16 and 1.35 10~ 16 W/m 2 for the 63 
and 145 /mi line, respectively. The upper limit for the 63 /mi 
line is only a factor 1.6 smaller than the flux predicted from the 
10~ 3 M disk model. The 145 //m upper limit is less sensitive and 
hence even consistent with our 10 -2 M disk model. Changing 
the outer disk radius to 500 AU does not affect the [0 1] lines. 

The observed [C 11] line is much stronger than in any of our 
models, suggesting indeed significant contamination from dif- 
fuse background material in the large ISO beam (80"). 
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Table 6. Observed fluxes for MWC480 



Type 


A 


Flux 


Reference 




63 jt/m 


< 1 1 HV 13 W/m-//im 
5s 1.1 1U W/1I1 /pill 


U ) 


Oi 


145 yum 


< 1.8 1CT 16 W/m 2 /pm 


(1) 


Cn 


158 pm 


5.9 ± 0.6 10- 16 W/m 2 


(1) 






5.5 ± 0.4 lO" 16 W/m 2 


(2) 


CO 


345 GHz 


2.88 K km/s 


(3) 






4.23 10~ 19 W/m 2 


(*) 


cont. 


1.3 mm 


360 ± 24 mjy 


(4) 


cont. 


2.7 mm 


39.9 ± 0.8 mjy 


(4) 



Refere nces : (1) | Cree ch-Eakma n~et al.l d2002[). (2) lLorenzetti et alj 
d2002h , (3) lThi et alJd2004h , (4) lMannings & Sargentld 19971) . 

(*) / T mb dv = W- 2 A c J(2k) (QJ- 1 / F,,(W/m 2 /Hz) dv, with Q a being 
the solid angle of the beam, i.e. 7r(13.7"/2) 2 for the JCMT. 

Our 2.2 10 -2 M model overpredicts the CO J=3-2 line flux 
by approximately one order of magnitude. By reducing the max- 
imum grain size, the total grain surface area increases and dust 
and gas temperatures couple in the region where this CO line 
forms; thus yielding (T g ) ~ 40 K and a factor two lower line 
fluxes. Another factor two can each be gained by reducing the to- 
tal disk mass and the disk outer radius to 500 instead of 700 AU. 

The observed continuum flux at 1 .3 mm is a factor 4 smaller 
than computed from the 2.2 10" 2 M disk model. Reducing 
again either the total disk mass or the maximum grain size, re- 
duces the modeled continuum fluxes to the right value. 

To conclude, by tuning the input parameters, the models can 
get close to the observed continuum as well as line fluxes, thus 
providing additional confidence that the models capture the es- 
sential physics and chemistry of protoplanetary disks. Hence, be- 
sides studying physics and chemistry of protoplanetary disks in 
a more general context, these models can also be used in the de- 
tailed analysis of individual objects that possess a larger set of 
continuum and line observations. 

7. Conclusions 

We used here a limited series of Herbig Ae disk models com- 
puted with the new thermo-chemical disk code called ProDiMo 
to study the origin and diagnostic value of the gas line tracers 
[Cn], [Oi] and CO. We do not include in this study effects of 
X-ray irradiation, grain settling or mixing. Even though X-rays 
are generally of minor importance for Herbig stars, grain set- 
tling and large scale mixing processes could affect the conclu- 
sions. The effect of grain settling will be included in a forthcom- 
ing larger model grid, but mixing processes need to be addressed 
with dynamical time-dependant models. The Monte-Carlo radia- 
tive transfer code Ratran is used to compute line profiles and 
integrated emission from various gas lines. 
The main results are: 

- The [C n] line originates in the disk surface layer where gas 
and dust temperatures are decoupled. The total line strength 
is dominated by emission from the disk outer radius. Thus 
the 157.7 fj.m line probes mainly the disk extension and outer 
disk gas temperature. The line forms in LTE. 

- The [O i] lines originate also in the disk surface layer even 
though somewhat deeper than the [Cn] line. The main con- 
tribution comes from radii between 30 and 100 AU. The [O i] 



lines form partially under NLTE conditions. Differences in 
line emission from escape probability and Monte Carlo tech- 
niques are smaller than 10%. 

- CO submm lines are optically thick down to very low disk 
masses of < 10 -4 M and form mostly in LTE. r gas = Jjust 
is not a valid approximation for these lines. Differences in 
line emission from escape probability and Monte Carlo tech- 
niques are smaller than 3%, except in the case of very opti- 
cally thin disk models (10~ 5 and 10~ 6 M ). 

- The [Oi] 63/145 pm and [Oi]63/[Cii] 158 //m line ratios 
trace disk mass in the regime between 10~ 2 and 10~ 6 M Q . 
Since the [C n] 158 fim line is very sensitive to the outer disk 
radius, the [O i] 63/[C n] 158 /mi is degenerate in that respect 
and its use requires additional constraints from ancilliary gas 
and/or dust observations. The sensitivity of these two line ra- 
tios to the dust grain sizes underlines the importance of using 
SED constraints along with the gas modeling to mitigate the 
uncertainty of dust properties. 

- A combination of the [Oi] 63/145 /zm and [Oi]63/[Cii] 158 
pm line ratios can be used to diminish the degeneracy caused 
by an unknown outer disk radius. 

- Neither total CO submm line fluxes nor line ratios can be 
used to measure the disk mass. However, the low rota- 
tional lines studied here provide an excellent tool to measure 
the disk outer radius and can thus help to mitigate the de- 
generacy between gas mass and outer radius found for the 
[Oi] 63/[Cn] 158 fim line ratio. 

Acknowledgements. We thank Dieter Poelman for fruitful discussions about ra- 
diative transfer methods and detailed code comparisons. 
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